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1 Abstract 

A meso-scale analysis is performed to determine the fracture process zone of concrete subjected to uniaxial 
tension. The meso-structure of concrete is idealised as stiff aggregates embedded in a soft matrix and 
separated by weak interfaces. The mechanical response of the matrix, the inclusions and the interface 
between the matrix and the inclusions is modelled by a discrete lattice approach. The inelastic response 
of the lattice elements is described by a damage approach, which corresponds to a continuous reduction 
of the stiffness of the springs. The fracture process in uniaxial tension is approximated by an analysis of 
a two-dimensional cell with periodic boundary conditions. The spatial distribution of dissipated energy 
density at the meso-scale of concrete is determined. The size and shape of the deterministic FPZ is 
obtained as the average of random meso-scale analyses. Additionally, periodicity of the discretisation is 
prescribed to avoid influences of the boundaries of the periodic cell on fracture patterns. The results of 
these analyses are then used to calibrate an integral-type nonlocal model. 



2 Introduction 

The energy dissipation during the fracture of concrete is influenced by the meso-structure of the material, 
the loading applied and the geometry of the specimen. The spatial distribution of the dissipated energy 
density across the fracture process zone (FPZ) is governed by both statistics and mechanics. Depending 
on the type of loading, either statistical or mechanical processes dominate. For instance, the shape of 
shear bands in concrete is strongly influenced by the mechanical interaction of aggregates. For this type 
of loading, the deterministic contribution on the shape of the FPZ is significant [30j. For tensile fracture 
the tortuosity of the path of the main crack, which dissipates most of the energy, is predominantly deter- 
mined by the statistics of the random arrangement of aggregates. Therefore, the crack paths in concrete 
subjected to the same loading conditions differ significantly, which implies that a purely deterministic 
meso-scale model, which does not consider the statistical variation of the fracture paths, cannot describe 
the FPZ of concrete subjected to tension. Hence, a direct determination of the mean FPZ by meso-scale 
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analysis requires averaging of the results of meso-scale analyses. Previous modelling efforts have been 
focused on an indirect determination of the FPZ by inverse analysis [1, 0> EI, 21, 24]. In one of these 



approaches, strength and fracture energies obtained from specimens of different sizes have been used to 
determine the width of the FPZ 27]. However, the results obtained with this calibration approach were 
not unique and provided only limited insight into the shape and size of the fracture process of concrete. 

The present study aims at determining the FPZ for mode-I fracture of concrete by meso-scale analysis. 
The nonlinear finite element method was used to analyse the spatial distribution of dissipated energy 
density on the meso-scale of concrete. The size and shape of the deterministic FPZ was determined as the 
average of meso-scale analyses. The results were then used to calibrate a deterministic nonlocal model, 
which can be used for large-scale structural analysis. To the authors' knowledge, this type of meso-scale 
analysis for the determination of the fracture process zone has not been carried out before. 

Within the framework of the nonlinear finite element analysis, three main approaches to fracture modeling 
can be distinguished. Continuum approaches describe the fracture process by higher-order constitutive 
models, such as integral-type nonlocal models [l|, [lj|. In continuum models with discontinuities, cracks 
are described as displacement discontinuities, which are embedded into the continuum description [26j |. 
Finally, discrete ap pro aches describe the nonlinear fracture process as failure of discrete elements, such as 
trusses and beams 29|, 11 1. In recent years, one type of discrete approach based on a lattice determined 



by Voronoi tesselation has been shown to be suitable for fracture simulations [f|. The lattice approach 
is robust, computationally efficient, and allows for fracture description by a stress- inelastic displacement 
relationship, similar to continuum models with discontinuities. With a specially designed constitutive 
model, the results obtained with this approach were shown to be mesh-independent |6[. Such a lattice 
approach is used for the meso-scale modelling in the present study. 

In the present work, discrete approaches are divided in the group of particle and lattice models. In particle 
models, the arrangement of particles can evolve, so that neighbours of particles might change during 
analysis. Therefore, particle models are suitable to describe processes involving large displacements. On 
the other hand, in lattice models the connectivity between nodes is not changed during the analysis, so 
that contact determination is not required. Consequently, lattice models are mainly suitable for analysis 
involving small strains 2(J 3, HI • 



In lattice analyses, the meso-structure of concrete can be described in at least two ways. In the first 
approach, the interaction between aggregates is modelled directly by single lattice elements [43| . All the 
nonlinearity of the material response between the aggregates is represented by the stress-strain response 
at a single point. This approach is characterised by computational efficiency, since the nodes of the 
finite element mesh correspond to the centres of concrete aggregates [ijj]- In the second approach, 
information on the heterogeneous meso-structure of concrete is mapped on a lattice in the form of spatially 
varying material properties [40j | . This approach requires a finer resolution, since individual aggregates are 
represented by several lattice elements. In the present study the latter approach is chosen, since a detailed 
description of the tortuous crack patterns is of importance. Because of the fine lattice used, the present 
study is limited to two-dimensional plane stress analyses with aggregates idealised as cylindrical inclusions. 
This oversimplifies the meso-structure of concrete so that a direct comparison with experimental results 
is difficult. However, the study may produce qualitative results that are useful for the development of 
macroscopic nonlocal constitutive models. 

The lattices used were generated from randomly placed nodes, which reduced the influence of the align- 
ment on the fracture patterns 23j, 39]. Similar observations have been made for other fracture approaches, 



in which irregular meshes can reduce the influence of the discretisation on fracture patterns [251 ] . The num- 
ber of degrees of freedom of the background lattice is reduced by placing lattice elements perpendicular 
to the boundary of aggregates Small aggregates, which cannot be captured by the background lattice, 
are approximated by a random field, which is generated by a spectral representation [42I fill]. Thus, the 
influence of the background lattice on the fracture patterns in the regions between the discretely modelled 
aggregates is reduced. 
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The constitutive response of the meso-structure of concrete can either be described by micro-mechanical 
models based on multi-scale analysis Q or, alternatively, by p henomenological constitutive models, which 
are commonly based on the theory of damage mechanics [311 ] , plasticity [13J and combinations of damage 
and plasticity [2^|. Plasticity and combinations of damage and plasticity are well suited to describe the 
fracture process of concrete in compression n the other hand, for mode-I fracture, damage 

constitutive models often provide satisfactory results [25| . In the present study an isotropic damage 
constitutive model is used. 



3 Meso-scale modelling approach 

The present approach to modelling the fracture process zone of concrete is based on a meso-scale de- 
scription. Aggregates, interfacial transition zones (ITZ) and mortar are modelled as separate phases with 
different material properties. For the mortar and ITZ, a random held of strength and fracture energy is 
applied. A lattice approach Q is used in combination with a damage mechanics model. 



3.1 Lattice approach 

The material response is represented by a discrete system of structural elements. The nodes of the lattice 
are randomly located in the domain, subject to the constraint of a minimum distance d m i n [43j |. For the 
discretely modelled aggregates, the lattice nodes are not placed randomly but at special locations, such 
that the middle cross-sections of the lattice elements form the boundaries between aggregates and mortar 
Q • The lattice elements are obtained from the edges of the triangles of the Delaunay triangulation of the 
domain (Fig. Hk), whereby the middle cross-sections of the lattice elements are the edges of the polygons 
of the dual Voronoi tesselation Q ■ 

Each lattice node possesses three degrees of freedom (two translations and one rotation). In the global 
coordinate system, shown in Figure [TJd, the degrees of freedom u c = (ui, V\, <fii,U2,V2, <f>2 ) T of the lattice 
nodes are linked to the displacement discontinuities u c = (u c ,v c ) T in the local co-ordinate system at 
point C, which is located at the center of the middle cross-section of the element. The relation between 
the degrees of freedom and the displacement discontinuities at C is 



Bu n 



(1) 



where 



B 



-cos a -sma — e cos a sin a e 
sin a —cos a —h/2 cos a sin a —h/2 



(2) 



The displacement discontinuities u c at point C are transformed into strains e = (e n , e s ) T = u c //i, where 
h is the distance between the two lattice nodes. The strains are related to the stresses <x = (cr n , a s ) T by 
an isotropic damage model, which is described below. 

The stiffness matrix of the lattice element in the local coordinate system has the form 

K = ^B T DB (3) 

where / is the cross-sectional area (for a two-dimensional model this area reduces to the length of the side 
shared by neighboring Voronoi polygons) and D is the material stiffness matrix. 

Heterogeneous materials are characterised by spatially varying material properties. In the present work 
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Figure 1: (a) Lattice based on Voronoi polygons, (b) Lattice element in the global coordinate system. 



this is reflected at two levels. Aggregates with diameters greater than </> m i n are modelled directly. The 
random distribution of the aggregate diameters (f> is defined by the cumulative distribution function used 
in 0, [13]. The aggregates are placed randomly within the area of the specimen, avoiding overlap of 
aggregates. Overlaps with boundaries are permitted. The heterogeneity represented by finer particles is 
described by an autocorrelated random field of tensile strength and fracture energy, which are assumed 
to be fully correlated. The random field is generated using a spectral representation [42|, El[ ■ This mixed 
approach, in the form of a discrete representation of the meso-structure and random field, is a compromise 
between model detail and computational time. 

The random field is characterised by an exponential autocorrelation function 

i?(0=exp(H£|7& 2 ) (4) 

and a Gaussian probability distribution function with a threshold to exclude negative values of strength 
and fracture energy. Parameter £ is the separation distance and parameter b is related to the autocorre- 
lation length l a as 

2/ 

(5) 



The autocorrelation length l a is independent of the spacing G? ln ; n of the lattice nodes [14 1. 
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3.2 Constitutive model 



An isotropic damage model is used to describe the constitutive response of ITZ and mortar. In the 
following section, the main equations of the constitutive model are presented. The stress-strain law reads 



cr = (1 — u>) D e = (1 — uj) cr 



(6) 



where cr = (cr n , a s ) T is the nominal stress, u> is the damage variable, D e is the elastic stiffness, e — (e n , e s ) T 
is the strain, and er = (<7 n ,CT s ) T is the effective stress. 



The elastic stiffness 



D, 



E 
jE 



(J) 



depends on model parameters E and 7, which control Young's modulus and Poisson's ratio of the material 
For a regular lattice and plane stress, Poisson's ratio is 



and the elastic stiffness is 



E u 



2E 



1-7 
3 + 7 

1 + 7 

3 + 7 



(8) 
(9) 



For a positive shear stiffness, i.e. 7 > 0, Poisson's ratio is limited to v < — , which is acceptable for 

o 

concrete but may be unrealistic for certain other materials. The damage parameter u> is a function of a 
history variable k, which is determined by the loading function 



/(e, k) = £ oq (e) - t 



(10) 



and the loading-unloading conditions 



/ < 0, k > 0, kf = 



(11) 



The equivalent strain 



eeq(en,e,) = ieo(l-c) + yQeo(c-l)+e n J + (12) 
corresponds to an elliptic stress envelope in the nominal stress space (Fig. [5]). 

For pure tensile loading, the nominal stress is limited by the tensile strength / t = Esq. For pure shear 
and pure compression, the nominal stress is limited by the shear strength / q = qf t and the compressive 
strength f c = eft, respectively. 

The softening curve of the stress-strain response in uniaxial tension is defined by the relation 



CT=/tex H"iJ (13) 

where w c = cuhe is considered as an equivalent crack opening under monotonic uniaxial tension. The 
stress-strain law in uniaxial tension can also be expressed in terms of the damage variable as 

<t = (1-lu)Ee (14) 
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Figure 2: Elliptic stress envelope in the nominal stress space. 



Comparing the right-hand sides of (|13p and (|14|1 . and replacing e by k, we obtain the equation 

( tj->hn\ 

(1 — uj) k = £o exp (15) 

V Wf / 

from which the damage parameter w can be determined iteratively. Parameter u>f determines the initial 
slope of the softening curve and is related to the meso- level fracture energy Gft = ftW{. The present 
damage constitutive model used here is conceptually similar to [13]. However, a more complex equivalent 
strain formulation involving the additional parameter c was required in the present study, since the 
compressive strength of concrete is significantly greater than its tensile strength. For c = 1 the original 
formulation proposed in [331 ] is obtained. 



3.3 Periodic cell 



One important issue in micro- and meso- mechanical modeling is the specific choice of boundary conditions 
at the boundary of the computational cell. In elastic homogenization theory, it is well known that 
kinematic boundary conditions (imposed displacements) lead to an apparent stiffness (describing the 
relation between average strain and average stress) that overestimates the effective one, while static 
boundary conditions (imposed tractions) lead to an underestimated stiffness. The influence of boundary 
conditions on the apparent elastic properties fades away with increasing size of the cell. If the cell is large 
enough, it is considered as a representative volume clement and the apparent clastic properties arc close 
to the effective properties of the homogenized medium. 

If the analysis is extended to highly nonlinear material behavior that leads to localization of inelastic 
processes, the boundary conditions may have a strong influence even if the cell is very large. The reason is 
that the imposed constraints usually lead to stress concentrations in a material layer near the boundary, 
which then bias the localization pattern. It is therefore preferable to use boundary conditions that are free 
of such bias and provide a statistically homogeneous distribution of localization patterns. The undesired 
boundary effects can be eliminated by the periodicity requirement that replaces the boundary conditions. 
This requires a special meso-structure and its discretization (in our case lattice), both compatible with 
the periodicity requirement. 

The meso-scale approach used in this paper is based on the concept of an initially rectangular computa- 
tional cell that can be periodically repeated in the directions parallel to its sides, both in the initial and 
the current (deformed) configuration. In addition to the displacements and rotations of individual lattice 
nodes that are contained withing the cell, the average strain components are considered as degrees of 
freedom of the computational model. Force quantities work-conjugate to these global degrees of freedom 



G 



Figure 3: Periodic cell. 



are directly related to the average stress components. It is thus easy to impose various types of mixed 
loading conditions. For instance, for uniaxial tension at the macroscopic scale, the average strain in the 
direction of loading is incremented (playing the role of the control variable) while the average stress in 
the lateral direction and the average shear stress are set to zero. 

When setting up the discretised equilibrium equations of the meso-scale model, the internal lattice ele- 
ments, i.e. those connecting two nodes that are both inside the cell, are processed in the standard way. 
Special treatment is needed for the boundary elements that connect one node inside the cell with another 
node physically located in one of the neighboring cells. The external lattice node is a periodic image of 
one of the internal nodes. For easier reference, let us denote I and J the internal lattice nodes and I' 
and J' their periodic images outside the cell; see Fig. [31 

The corresponding lattice element can be considered as connecting either nodes I and J', or I' and J. 
Both representations arc equivalent and only one of them, say I and J', is actually processed when setting 
up the contribution to the equilibrium equations. Due to the assumption of periodicity, the rotations of 
lattice nodes J and J' are the same and the translations of node J' can be expressed as the translations 
of node J plus an appropriate contribution of the average strains. More specifically, suppose that the 
initial position of node J' is shifted with respect to node J by k x a in the x-direction and by k y b in the 
y-direction, where a and b are dimensions of the cell and k x and k y are integers —1, or 1 (for the specific 
case plotted in Fig. OH we have k x = 1 and k y = 0). The translations of node J' are then expressed as 

Ujt = uj + k x aE x + k y bE xy (16) 
vj> = vj + kybE y (17) 

where E x and E y are the average normal strains in the x- and y-direction, respectively, and E xy is the 
average engineering shear strain. Note that the contribution of the average shear has been attributed 
exclusively to the translations along the x-direction. This is perfectly justified, because rigid-body rotation 
of the entire lattice is completely arbitrary and we can fix it by assuming that straight lines parallel to 
the x-axis do not rotate. The other two rigid-body modes (translations) need to be suppressed by setting 
the translations of one selected lattice node to zero. 

Making use of <P^ j) — (jT7)) and of the relation for the rotations, <f>j/ = <j>j, we can set up the transformation 
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Table 1: Model parameters. 



E [GPa] 
Matrix 30 
ITZ 45 
Aggregate 90 



7 ft [MPa] 
0.33 5.3 
0.33 1.8 
0.33 



q c G ft [N/m 
2 10 93 
2 10 31 
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When this rule is combined with equation (JXJ) linking the displacements of nodes / and J' and the 
discontinuities, the 6x9 transformation matrix from (| 1S|) . denoted as T, multiplies matrix B from the 
right. It follows from duality that the internal forces must be multiplied by T T from the left before 
they are inserted into the equilibrium conditions. The usual 6x6 stiffness matrix of the element, K, 
is therefore transformed into a 9 x 9 matrix T T KT before it enters the assembly process. The global 
stiffness matrix remains banded, with the exception of the three rows and columns that contain terms 
related to the global degrees of freedom (average strains) and to the corresponding equilibrium equations 
(related to the average stresses). 



4 Meso-scale analysis 

4.1 Geometry, loading setup and results 

The aim of the present study is to analyse the fracture process zone of concrete under uniaxial tension. 
A specimen with periodic boundary conditions, material properties and background lattice was used to 
reduce the influence of boundaries on the fracture process. The material properties for the mortar, ITZ 
and aggregate phases were chosen according to Tab.[TJ Parameters ft and G{ in Tab.[T]are the mean values 
of the random fields of tensile strength / t and fracture energy Gf , respectively. The material parameters 
were chosen to approximate a macroscopic tensile strength of 3 MPa and a macroscopic fracture energy 
of 100 J/m 2 . The ratio of the Young's modulus of matrix and aggregate is assumed to be 3. The Young's 
modulus of the ITZ is 

E lT z = 1 2 1 (19) 
E m E a 

where E m and E a are the Young's moduli of matrix and aggregate, respectively. Furthermore, the ratios 
of tensile strength and fracture energy for matrix and ITZ are 3. These values are in the range used for 



other meso-scale analyses in [17 



The specimen used in the analyses (Fig. 2]) had a height and width of a = b = 100 mm, respectively. 
The lattice for the meso-scale analysis was generated with d m i n = 0.75 mm. The aggregate volume 
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Figure 4: (a) Geometry and loading setup of the periodic specimen with a = b = 100 mm. (b) Load- 
displacement curve for one random meso-scale analysis. 



fraction was chosen as p = 0.3 with a maximum and minimum aggregate diameter of <fi max = 12 mm 
and (/> m in = 4.75 mm, respectively. The approach to generate the distribution of aggregate diameters 
is described in [l7j |. Furthermore, the random field was characterised by the autocorrelation length 
Z a = 1 mm and the coefficient of variation c v = 0.2. The specimen was subjected to uniaxial tensile 
stress <7y, which was controlled by the displacement in y-direction. On the left and right boundary of the 
specimen, boundary conditions that correspond to a vanishing macroscopic stress in the ^-direction were 
applied. 

The stress-strain response for one random meso-scale analysis is shown in Fig. [4]b. The crack patterns for 
three stages of the uniaxial tensile loading, marked in Fig. 0b, are shown in Fig. [51 Red (dark grey) lines 
mark cross-sections of elements in which the damage parameter increases at this stage of analysis (i.e., 
active, opening cracks), whereas light grey lines indicate cross-sections in which the maximum damage 
parameter was reached at an earlier stage of analysis (i.e., passive, closing cracks). From Fig. Eh and it 
can be seen that the path of the active crack is highly tortuous. At earlier stages of analysis, distributed 
cracking occurs, whereby the cracks are located mainly in the ITZs (light grey cracks) (Fig. [5] a). Shortly 
after peak (Fig. 0b), many of these cracks cease to grow, i.e. their opening decreases, and only few cross- 
sections exhibit increasing crack opening (red or dark grey lines) . This is in agreement with observations 
made in Energy dissipation in these localised cracks constitutes a substantial part of 

dissipated energy. In Fig. [6j the mean of 100 meso-scale analyses is presented with error bars showing 
the standard deviation. The results are plotted in terms of the average stress (force per unit length of 
the unit cell side and per unit thickness) against the average strain (relative displacement of the opposite 
sides divided by their distance). The average stress-strain curve exhibits pre-peak nonlinearities typical 
of uniaxial tensile experiments of concrete. The standard deviation is small in the pre-peak regime, but 
increases strongly in the post-peak regime of the analysis. 

In addition to the stress-strain curves in Fig. the dissipated energy densities are evaluated. The 
averaging of dissipated energy densities is complicated, since the location of the final crack, in which most 
of the energy is dissipated, depends on the meso-structure of the material and cannot be determined in 
advance. Direct averaging of the energy densities for random meso-scale analyses, would lead, in the limit 
of an infinite number of analyses, to a uniformly distributed energy density over the specimen height. 
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Figure 5: Crack patterns for three stages of loading marked in Fig. UJa. Red (dark grey) lines indicate 
cross-sections with increasing crack opening. Light grey lines indicate cross-sections with decreasing crack 
openings. 
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Figure 6: Mean curve obtained by averaging the results of 100 analyses. Error bars show the range 
between the mean plus and minus one standard deviation. 
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Figure 7: Mean profile of the dissipated energy density across the process zone obtained by averaging 
100 analyses with a resolution of 1.6 mm. Error bars show the range between the mean plus and minus 
one standard deviation. Fracture process zones are shifted before averaging. 

For a meaningful characterization of the random fracture process zones, the results of individual analysis 
were post-processed. Firstly the y-coordinate of the centre of the fracture process was determined by 
processing the dissipated energy of all elements. Then, all elements were shifted in the y-direction 
(Fig. U]), such that the y-coordinate of the centre of the dissipated energy density of each individual 
analysis coincides with the y-coordinate of the centre of the specimen. This shift is admissible, since the 
specimen is fully periodic, i.e. not only the boundary conditions, but also meso-structure and background 
lattice are periodic. In the next step, the specimen was subdivided in a regular rectangular grid of 
cells with x- and y-edge lengths of c x = a/64 and c y = &/64. The energy densities in these cells were 
determined by integrating the dissipated energy of all elements located within them, and by dividing these 
energies by the cells areas. Subsequent averaging of the energy densities of the 100 analyses results in 
the average dissipated energy density, which characterizes the fracture process zone. The energy density 
along the y-direction, i.e. perpendicular to the crack, is presented in Fig. [7] by averaging the dissipated 
energy along the crack (along the cc-direction) . The error bars in Fig. [7] represent the range between 
the mean plus and minus one standard deviation. The fracture process zone has its maximum value 
in the centre with symmetrically decreasing slopes on both sides. The width of the fracture process 
zone is mainly determined by the tortuosity of the cracks and is for the present meso-structure equal 
to approximately. 3 times the size of the maximum aggregates. The dissipated energy density in the 
x-direction (along the crack) is shown in Fig. [8l For this representation, the density is integrated in the 
y-direction. Thus, the result corresponds to the energy dissipated per unit area of the ligament (per unit 
length in the two-dimensional setting). The distribution of dissipated energy along the crack in Fig. [5] 
is fluctuating around a constant mean value, which corresponds to the fracture energy of the material. 
Again, error bars indicate the mean plus and minus one standard deviation. 

4.2 Discussion of the meso-scale analysis results 

The standard deviations of the dissipated energy densities are considerably greater than those of the 
stress-strain curves. At some points of the fracture process zone in Fig. [7| the standard deviation is even 
greater than the mean value, which indicates that the statistical distribution is far from normal (since all 
dissipation values are non-negative). The difference in standard deviation is expected, since the average 
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Figure 8: Mean values of the dissipated energy density along the process zone obtained by averaging 100 
analyses with a resolution of 1.6 mm. Error bars show the range between the mean plus and minus one 
standard deviation. 



stress-strain curves are global results, which involve averaging of local responses, whereas the dissipated 
energy density is local. With such significant standard deviations, it is necessary to check whether 100 
analyses are sufficient to obtain statistically representative results. Therefore, two additional sets of 100 
analyses were carried out. The averages of the three sets are presented in Fig. [§1 and [TU1 in the form of 
stress-strain curves and fracture process zones, respectively. The difference between the results of the 
three sets is small, which indicates that 100 analyses are sufficient to obtain statistically representative 
results. 

Additionally, the influence of the lattice on the results was investigated. In the present study, the size of 
lattice elements and their spatial arrangement was chosen independently of the material meso-structure. 
Of course, the lattice elements must be sufficiently small to be able to represent the discretely modelled 
aggregates. However, fracture is modelled as displacement jumps within individual lattice elements. 
Therefore, the spatial orientation and size of the lattice elements might influence the fracture patterns. 
The influence of these two parameters is investigated by two additional studies. Firstly, analyses with 
three different random lattices with the same minimum distance d m i n were conducted. The same aggregate 
arrangement and random field was used. The stress-strain curves and crack patterns of these analyses 
are compared in Fig. [TTIandfT^l The stress-strain curves exhibit only small differences in the post-peak 
regime, which are smaller than the standard deviations of the stress-strain curves shown in Fig. [51 The 
overall crack patterns in Fig. 1121 which determine the shape of the fracture process zone in Fig. [7J are 
almost independent of the background lattice. Only small differences are visible, which have a negligible 
influence on the tortuosity of the fracture process zones obtained. 

As the next step, the influence of the lattice element size was studied. In the constitutive model used in the 
present study, fracture is described by a stress-crack opening curve. Therefore, the results are expected 
to be independent of the size of the lattice elements. Three analyses with the same meso-structure 
but different background lattices with minimum distances d m i n = 1, 0.75 and 0.5 mm, respectively, 
were performed. The differences among the stress-strain curves (Fig. I13[) are small, which shows that 
the present lattice approach is free of pathological mesh dependence. Also, the crack patterns in Fig. [T3] 
exhibit only local differences, which are due to the random arrangement of lattice elements (compare with 
Fig. fT2|) . but not their size. These local differences are not expected to influence the average fracture 
process zones presented in Fig. [7J The two studies above show that discretely modelled aggregates 
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Figure 9: Comparison of average stress-strain curves obtained from three different sets of 100 analyses. 
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Figure 10: Comparison of average fracture process zones obtained from three different sets of 100 analyses. 
Fracture processes zones are shifted before averaging. 
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Figure 11: Comparison of average stress-strain curves obtained using three different discretisations with 
the same meso-structure and lattice element size, but different random positions of lattice nodes. 




Figure 12: Comparison of crack patterns obtained from analyses using three different lattices with the 
same meso-structure and the same average lattice element size. 
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Figure 13: Comparison of average stress-strain curves obtained from lattices with minimum distances 
of e? m i n = 1, 0.75 and 0.5 mm, respectively. In all three analyses the same aggregate arrangement and 
random field is used. 




Figure 14: Comparison of crack patterns obtained from lattices with minimum distances of d m in — 1 mm, 
rfmin = 0.75 mm and d m i n = 0.5 mm, respectively. In all three analyses, the same aggregate arrangement 
and random field is used. 

in connection with a random field for the mortar provide results which are almost independent of the 
background lattice. 

A periodic cell was chosen in the present study so that the results are independent of the boundary 
conditions and the fracture process zones can be shifted before the averaging of the results. The size of 
the periodic cell was chosen so that the tortuosity of the fracture patterns is statistically representative. A 
decrease of the size of the periodic cell in the direction of loading is expected to influence the tortuosity of 
the crack patterns, i.e. the width of the fracture process zone. The smaller the cell size, the less tortuous 
are the crack patterns. However, it is expected that there is an upper limit of the cell size, above which 
the width of the fracture process zone remains almost constant. The influence of the cell size was studied 
by performing two sets of 100 analysis for two additional cells with a = 50 and 150 mm. The average 
stress-strain curves and fracture process zones for these two sizes are compared to the original results for 
a = 100 mm in Figs. fTSI and I16[ respectively. Although the difference between the average stress-strain 
curves for the three sizes is small, it can be seen that a reduction of the size leads to a steeper softening 
branch in the post-peak regime of the stress-strain curves. Furthermore, the comparison of the fracture 
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Figure 15: Comparison of average stress-strain curves from 100 analyses obtained for periodic cells of 
size (perpendicular to the loading direction) a = 50, 100 and 150 mm. 
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Figure 16: Comparison of average fracture process zones obtained for periodic cells of sizes a — 50, 100 
and 150 mm. Fracture process zones are shifted before averaging. 
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process zones shows that the width of the fracture process zone reduces considerably when the size of 
the periodic cell is decreased from 100 to 50 mm. On the other hand, an increase of the size from 100 to 
150 mm has only a small influence on the width of the fracture process zone. Consequently, a cell size of 
100 mm was deemed to be sufficient to represent the fracture process zone. 



5 Comparison of meso-scale analyses to macroscopic nonlocal 
constitutive model 

In the previous two sections, the meso-scale description of fracture in concrete was used to determine the 
average of dissipated energy densities. The present study adopts the view that this average of densities 
corresponds to the fracture process zone, which is represented by macroscopic nonlocal models. In the 
following section, a one-dimensional macroscopic nonlocal damage model is compared to the results of 
meso-scale analyses. 



5.1 Nonlocal damage model 

In the one-dimensional setting, the stress-strain law used by the damage model is 

a = (l-Lu)Ee (20) 

where a is the uniaxial stress, to is the damage variable, E is Young's modulus and e is the strain. Damage 
evolution is driven by the nonlocal strain s, which represents a weighted spatial average of the local strain 
e (in multiple dimensions, a scalar measure of the strain level called the equivalent strain would be used). 
In an unbounded one-dimensional medium, the nonlocal strain is evaluated as 

/oo 
a(x-0e(0^ (21) 
-oo 

where a is the nonlocal weight function, decaying with increasing distance between points x and £ and 
normalized such that the averaging operator does not modify a uniform field. The weight function is 
sometimes taken as a Gauss-type function, but in the present study we consider the quartic polynomial 
with a bounded support, 

a ( r ) = J5_ A _ 4\ (22) 
w 16 i? \ R 2 1 K ' 

and the exponential function 

a(r)=i«p(-M) (23) 

Both weight functions are normalized such that J_ a(r) dr — 1. Parameters R or I reflect the char- 
acteristic length of the material and have a strong influence on the width of the process zone, as will 
be shown later. Function (|2"3"f corresponds to Green's function of the differential equation e — l 2 e" = e 
(on an infinite one-dimensional domain). With this s pec ific choice, the nonlocal integral-type model is 
equivalent to the implicit gradient model proposed in [361 ] - 

Since damage is irreversible, the damage variable lu is related to an internal variable k, which represents 
the maximum level of nonlocal strain ever reached in the previous history of the material. Formally, n 
can be defined by the loading-unloading conditions (jlip with the loading function 

f(e, K )=e- K (24) 
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Table 2: Parameters that control the post-peak part of the stress-strain curve. 



parameter 


weight 


length 


£l 


£2 


£3 


£f 


set 


function 


parameter 


[io- 3 ] 


[io- 3 ] 


[10~ 3 ] 


[10- 3 ] 


A 




R = 20 mm 


0.255 


5 


0.0919 


0.603 


B 


mi 


i? = 15 mm 


0.238 


7 


0.0942 


0.814 


C 


422) 


i? = 10 mm 


0.221 


15 


0.0958 


1.335 


D 


(El 


^ = 8 mm 


0.280 


2.5 


0.0875 


0.451 


E 


(El 


I — 6 mm 


0.255 


3.5 


0.0919 


0.603 


F 


HMD 


I = 4 mm 


0.230 


7 


0.0950 


0.990 



The damage law that links the strain to damage is closely related to the shape of the stress-strain curve. 
Since the curves obtained from the lattice model exhibit nonlincarity already before the peak, we consider 
a damage law that can reproduce such behavior: 



1 — exp — 



1 



TO \ £ 



u = g{n) = < 



1 £ 3 

1 exp 



£r 



(^) 



if k < £1 
if k > £1 



(25) 



The primary model parameters are the uniaxial tensile strength / t , the strain at peak stress (under 
uniaxial tension) e p , and additional parameters e±, £2 and n, which control the post-peak part of the 
stress-strain law. Other parameters that appear in (EH} can be derived from the condition of zero slope 
of the stress-strain curve at k — £ p and from the conditions of stress and stiffness continuity at n = £\: 



1 



rn 



£r 



£3 



ln(££ p // t ) 
£1 

(£l/£p) m - 1 
f 1 

£1 exp — 

V m \ £ D 



(26) 
(27) 

(28) 



5.2 Parameter identification 

The parameters have been adjusted so as to fit the uniaxial stress-strain curve in Fig. [S]and the fracture 
process zone in Fig. [71 The basic parameters E — 29.6 GPa, / t = 2.86 MPa and £ p = 0.198 x 10~ 3 can be 
directly determined from the ascending branch of the stress-strain curve, and the corresponding value of 
m = 1.39 is calculated from (|26p . For the one-dimensional model, the strain profile remains uniform up 
to the peak, and so the above-mentioned parameters are independent of the characteristic length of the 
nonlocal model. The remaining parameters are related to the descending branch and their values must 
be optimized for each specific choice of the nonlocal weight function separately. Here we consider two 
different nonlocal weight functions, (|22[) and (ESI), each with three different values of the characteristic 
length parameter R or I. The resulting six sets of parameters are listed in Table El Parameter n controls 
the shape of the tail of the stress-strain curve and can be taken by the same value n = 0.85 for all the 
parameter sets considered here. Parameters E\ and £2 were obtained by fitting and parameters £f and £3 
were then calculated from ([2"Tj) - (l2"51) . 

The comparison of the nonlocal model and the meso-scale approach is shown in Figs. [T7] and [TBI With a 
proper choice of parameters of the macroscopic nonlocal model, the damage evolution law (|25|) allows for 
a good fit of the average stress-strain curve obtained from the meso-scale lattice analyses; see Fig.[T7] The 
agreement is almost perfect for all six sets of parameters presented in Table [2J This means that the same 
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Figure 17: Comparison of average stress-strain curves obtained with the meso-scale lattice model and 
with the macroscopic nonlocal model using different combinations of parameters. 

global response can be reproduced with different types of nonlocal weight functions and different values 
of the length parameter R or I. The optimal value of the characteristic length cannot be deduced from 
the average stress-strain curve, which is in fact the rescaled load-displacement curve and characterizes 
the global response. As shown in Fig. [TBI the characteristic length controls the width of the process zone. 

On the other hand, the specific form of the weight function has only a weak influence on the distribution 
of energy dissipation across the process zone; see Fig. 1191 The exponential weight function gives a slightly 
higher peak of the dissipation profile but this difference is much less pronounced than the difference in the 
weight functions; see Fig. [2UJ Of course, the length parameters R and I do not have the same meaning; 
similar dissipation profiles are obtained if R is set approximately to 2.5/. For R = 15 mm or I = 6 mm, 
a good fit of the dissipation profile predicted by the meso-mechanical lattice model is achieved, except 
for the center of the process zone where the extremely high dissipation density obtained with the lattice 
model is underestimated. 



6 Conclusions 

In the present work, a meso-scale approach was used to determine the fracture process zone of concrete, 
which is characterized by the average of dissipated energy densities. Then, a macroscopic nonlocal model 
was fitted to the meso-scale results. The work resulted in the following conclusions: 

• The average of dissipated energy densities obtained from the meso-scale analyses resulted in a 
fracture process zone of a finite width, which is determined by the tortuosity of the crack path. 

• The meso-scale results obtained with a lattice approach are independent of the size and spatial 
orientation of the elements. This is achieved by describing the material properties by a random 
field of strength. 

• The standard deviation of the energy densities obtained is greater than the mean value, which 
indicates a deviation from the normal distribution. Nevertheless, the results were shown to be 



lattice model 
nonlocal model, parameter set A 
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Figure 18: Distributions of the dissipated energy density across the fracture process zone obtained with 
the meso-scale lattice model and with the macroscopic nonlocal model using (a) quartic polynomial weight 
function (f22|) . (b) exponential (Green- type) weight function (|23|) 
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Figure 19: Distributions of the dissipated energy density across the fracture process zone obtained with 
the macroscopic nonlocal model using quartic polynomial weight function (|22[) and exponential (Green- 
type) weight function ([2"5|) 
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Figure 20: Quartic polynomial weight function (f22l) and exponential (Green-type) weight function (|23|) 
with the ratio of length parameters R/l — 2.5. 
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statistical representative. 



In future work, the meso-scale results will be used to calibrate nonlocal macroscopic damage models 
reported in the literature, which describe boundaries differently. These nonlocal models will then be 
applied to the analysis of notched concrete beams, for which boundary conditions play an important 
role. Additionally, meso-scale analysis of these concrete beams will be carried out, which will allow the 
identification of nonlocal models suitable for the description of the failure of concrete. 
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